Our primary hypotheis testing method is to fit generalized linear mixed model that includes the factor “generation” as a fixed effect and then determine if the estimated effect of the different levels of this factor are significantly different from one another, once all other effects have been controlled for. This is likely the best way to address our questions given that the covariates length and year have such large effects on fitness, but let’s not forget the clever approach that has been employed in past studies of RRS.
The delta method takes advantage of the fact while TLF is a count variable with variance greater than the mean and requiring very advanced statistics to appropriately model, RRS is a ratio. If we are able to describe mean and variance of TLF for each of the groups we should be able to divide them by one another and generate a ratio, and a confidence interval for this ratio can be estimated using the a bit of clever math and likelihoods. If the confidence interval doesn’t include one, then the fitness of the groups is different!
Let’s use this approach for each year.
Overall RRS
# kalinowski method
# rather than re-invent the wheel here let's make sure my interpretation of the Kalinowski paper is correct and modify seom existing code from a reviewed paper. My understanding matched the modified code here, so I just used that. We should be sure to cite if it winds up in a paper.
# https://doi.org/10.1098/rsos.221271
rrs_ci_kalinowski_auke <- function(n_h_off, n_w_off, n_h_par, n_w_par, alpha){
chi_alpha <- qchisq(p = (1 - alpha), df = 1)
n_off <- sum(c(n_h_off, n_w_off))
n_par <- sum(c(n_h_par, n_w_par))
rs_h <- n_h_off / n_h_par
rs_w <- n_w_off / n_w_par
p_h_par <- n_h_par / n_par
p_w_par <- n_w_par / n_par
rrs_h <- rs_h / rs_w
rrs_w <- rs_w / rs_w
rrs_avg <- (rrs_h * p_h_par) + (rrs_w * p_w_par)
rrs_ml <- (n_h_off * log(p_h_par * rrs_h / rrs_avg)) + (n_w_off * log(p_w_par * rrs_w / rrs_avg))
xi_dist <- bind_rows(
lapply(seq(0.01, 50, by = 0.01), function(rrs_h_xi) { #upper RRS up to 50
rrs_avg_xi <- (rrs_h_xi * p_h_par) + (rrs_w * p_w_par)
tibble(rrs_crit = rrs_h_xi,
logl = (n_h_off * log(p_h_par * rrs_h_xi / rrs_avg_xi)) + (n_w_off * log(p_w_par * rrs_w / rrs_avg_xi)) - (rrs_ml - chi_alpha / 2)
)
} )
)
rrs_min <- xi_dist %>%
mutate(abs_logl = abs(logl)) %>%
filter(rrs_crit < rrs_h) %>%
top_n(-1, abs_logl) %>%
pull(rrs_crit)
rrs_max <- xi_dist %>%
mutate(abs_logl = abs(logl)) %>%
filter(rrs_crit > rrs_h) %>%
top_n(-1, abs_logl) %>%
pull(rrs_crit)
return(c(rrs_min, rrs_h, rrs_max))
}
# get data to enter
rrs_delta <- F12_mmdata %>%
group_by(generation, year) %>%
summarise(n = n(), n_offspring = sum(tlf)) %>%
pivot_wider(names_from = generation, values_from = c(n, n_offspring))
#n_h_off, n_w_off, n_h_par, n_w_par, alpha
# embarrisingly redundant, slow code, consider revising, can be made at least 3x faster by not calling the same function 3 times witht e same parameters...
# argument order reminder of function, n_h_off, n_w_off, n_h_par, n_w_par, alpha
rrs_delta %<>%
rowwise() %>%
mutate(rrs_NH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[1],
rrs_NH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[2],
rrs_NH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[3]) %>%
mutate(rrs_FH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[1],
rrs_FH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[2],
rrs_FH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[3]) %>%
mutate(rrs_FN_lwr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[1],
rrs_FN = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[2],
rrs_FN_upr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[3])
Now let’s plot these.
plot_data <- rrs_delta %>%
select(year, starts_with("rrs")) %>%
pivot_longer(cols = -c(year), names_prefix = "rrs_", names_to = c("numerator", "value_type"), names_sep = "_", values_to = "value") %>%
mutate(value_type = case_when(is.na(value_type) ~ "MeanRRS",
TRUE ~ value_type)) %>%
pivot_wider(names_from = "value_type", values_from = "value") %>%
mutate(numerator = as.factor(numerator),
numerator = fct_relevel(numerator, "FN", "FH", "NH"))
ggplot(plot_data)+geom_errorbar(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, color = "darkgrey")+theme_bw()+geom_point(aes(x = year, y = MeanRRS, color = numerator), size = 3, position = position_dodge(width = 0.3))+xlab("")+scale_color_manual(name = "Contrast", values = c("#66CCEE", "#CCBB44", "#AA3377"), labels = c(expression(F[1]*" / NOR"),expression("HOR / " *F[1]), "HOR / NOR" ))+ylab(expression(" "[Delta]*"RRS"))+ theme(legend.text.align = 0)

ggplot(plot_data)+geom_errorbar(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, color = "darkgrey")+theme_bw()+geom_point(aes(x = year, y = MeanRRS, color = numerator), size = 3, position = position_dodge(width = 0.3))+xlab("")+scale_color_manual(name = "Contrast", values = c("#66CCEE", "#CCBB44", "#AA3377"), labels = c(expression(F[1]*" / NOR"),expression("HOR / " *F[1]), "HOR / NOR" ))+ylab(expression(" "[Delta]*"RRS"))+ theme(legend.text.align = 0) + coord_cartesian(ylim=c(0, 2))

#ggplot(plot_data)+geom_pointrange(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, color = "red")+theme_bw()
For the HOR RRSs, same results within each year (except 2015) using the delta method. RRS < 1 for regardless of whether F1s or NOR immigrants are in the numerator. Note that that this RRS calculation can’t take differences in covariates like length and release day into effect.
Also note one important difference F1/NOR fitness is less than one in 2012.
#Let's also calculate an overall RRS. Note that because some years have higher fitness than others, the unbalanced sample sizes across generation and year will create some major biases here and we shouldn't present this.
rrs_delta_all <- F12_mmdata %>%
group_by(generation) %>%
summarise(n = n(), n_offspring = sum(tlf)) %>%
pivot_wider(names_from = generation, values_from = c(n, n_offspring))
rrs_delta_all %<>%
rowwise() %>%
mutate(rrs_NH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[1],
rrs_NH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[2],
rrs_NH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[3]) %>%
mutate(rrs_FH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[1],
rrs_FH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[2],
rrs_FH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[3]) %>%
mutate(rrs_FN_lwr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[1],
rrs_FN = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[2],
rrs_FN_upr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[3])
Sex Specific RRS
We also chose not to include sex specific RRS because sex and sex*generation interaction were not retained in the final model, but let’s also split the delta method results up by sex and plot again to keep with conventions.
Males
# get data to enter
rrs_delta_males <- F12_mmdata %>%
filter(sex == "M") %>%
group_by(generation, year) %>%
summarise(n = n(), n_offspring = sum(tlf)) %>%
pivot_wider(names_from = generation, values_from = c(n, n_offspring))
## `summarise()` has grouped output by 'generation'. You can override using the
## `.groups` argument.
#n_h_off, n_w_off, n_h_par, n_w_par, alpha
# embarrisingly redundant, slow code, consider revising, can be made at least 3x faster by not calling the same function 3 times witht e same parameters...
# argument order reminder of function, n_h_off, n_w_off, n_h_par, n_w_par, alpha
rrs_delta_males %<>%
rowwise() %>%
mutate(rrs_NH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[1],
rrs_NH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[2],
rrs_NH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[3]) %>%
mutate(rrs_FH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[1],
rrs_FH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[2],
rrs_FH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[3]) %>%
mutate(rrs_FN_lwr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[1],
rrs_FN = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[2],
rrs_FN_upr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[3])
plot_data <- rrs_delta_males %>%
select(year, starts_with("rrs")) %>%
pivot_longer(cols = -c(year), names_prefix = "rrs_", names_to = c("numerator", "value_type"), names_sep = "_", values_to = "value") %>%
mutate(value_type = case_when(is.na(value_type) ~ "MeanRRS",
TRUE ~ value_type)) %>%
pivot_wider(names_from = "value_type", values_from = "value") %>%
mutate(numerator = as.factor(numerator),
numerator = fct_relevel(numerator, "FN", "FH", "NH"))
ggplot(plot_data)+geom_errorbar(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, color = "darkgrey")+theme_bw()+geom_point(aes(x = year, y = MeanRRS, color = numerator), size = 3, position = position_dodge(width = 0.3))+xlab("")+scale_color_manual(name = "Contrast", values = c("#66CCEE", "#CCBB44", "#AA3377"), labels = c(expression(F[1]*" / NOR"),expression("HOR / " *F[1]), "HOR / NOR" ))+ylab(expression(" "[Delta]*"RRS"))+ theme(legend.text.align = 0) +ggtitle("RRS Males")

#ggplot(plot_data)+geom_pointrange(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, col
#kinda difficult to see 2012 CI let's zoom in
ggplot(plot_data)+geom_errorbar(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, color = "darkgrey")+theme_bw()+geom_point(aes(x = year, y = MeanRRS, color = numerator), size = 3, position = position_dodge(width = 0.3))+xlab("")+scale_color_manual(name = "Contrast", values = c("#66CCEE", "#CCBB44", "#AA3377"), labels = c(expression(F[1]*" / NOR"),expression("HOR / " *F[1]), "HOR / NOR" ))+ylab(expression(" "[Delta]*"RRS"))+ theme(legend.text.align = 0) +ggtitle("RRS Males")+coord_cartesian(ylim=c(0,4))

Females
# get data to enter
rrs_delta_females <- F12_mmdata %>%
filter(sex == "F") %>%
group_by(generation, year) %>%
summarise(n = n(), n_offspring = sum(tlf)) %>%
pivot_wider(names_from = generation, values_from = c(n, n_offspring))
## `summarise()` has grouped output by 'generation'. You can override using the
## `.groups` argument.
#n_h_off, n_w_off, n_h_par, n_w_par, alpha
# embarrisingly redundant, slow code, consider revising, can be made at least 3x faster by not calling the same function 3 times witht e same parameters...
# argument order reminder of function, n_h_off, n_w_off, n_h_par, n_w_par, alpha
rrs_delta_females %<>%
rowwise() %>%
mutate(rrs_NH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[1],
rrs_NH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[2],
rrs_NH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_NORimmigrant, n_HOR, n_NORimmigrant, alpha = 0.05)[3]) %>%
mutate(rrs_FH_lwr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[1],
rrs_FH = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[2],
rrs_FH_upr = rrs_ci_kalinowski_auke(n_offspring_HOR, n_offspring_F1, n_HOR, n_F1, alpha = 0.05)[3]) %>%
mutate(rrs_FN_lwr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[1],
rrs_FN = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[2],
rrs_FN_upr = rrs_ci_kalinowski_auke(n_offspring_F1, n_offspring_NORimmigrant, n_F1, n_NORimmigrant, alpha = 0.05)[3])
plot_data <- rrs_delta_females %>%
select(year, starts_with("rrs")) %>%
pivot_longer(cols = -c(year), names_prefix = "rrs_", names_to = c("numerator", "value_type"), names_sep = "_", values_to = "value") %>%
mutate(value_type = case_when(is.na(value_type) ~ "MeanRRS",
TRUE ~ value_type)) %>%
pivot_wider(names_from = "value_type", values_from = "value") %>%
mutate(numerator = as.factor(numerator),
numerator = fct_relevel(numerator, "FN", "FH", "NH"))
ggplot(plot_data)+geom_errorbar(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, color = "darkgrey")+theme_bw()+geom_point(aes(x = year, y = MeanRRS, color = numerator), size = 3, position = position_dodge(width = 0.3))+xlab("")+scale_color_manual(name = "Contrast", values = c("#66CCEE", "#CCBB44", "#AA3377"), labels = c(expression(F[1]*" / NOR"),expression("HOR / " *F[1]), "HOR / NOR" ))+ylab(expression(" "[Delta]*"RRS"))+ theme(legend.text.align = 0) +ggtitle("RRS Females")

#ggplot(plot_data)+geom_pointrange(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, col
Hard to see given the very high range in 2015. We’ll make a second plot.
ggplot(plot_data)+geom_errorbar(aes(x = year, y = MeanRRS, ymin = lwr, ymax = upr, color = numerator), width = 0.3, position = position_dodge())+ geom_hline(aes(yintercept = 1), linetype = 2, color = "darkgrey")+theme_bw()+geom_point(aes(x = year, y = MeanRRS, color = numerator), size = 3, position = position_dodge(width = 0.3))+xlab("")+scale_color_manual(name = "Contrast", values = c("#66CCEE", "#CCBB44", "#AA3377"), labels = c(expression(F[1]*" / NOR"),expression("HOR / " *F[1]), "HOR / NOR" ))+ylab(expression(" "[Delta]*"RRS"))+ theme(legend.text.align = 0) +ggtitle("RRS Females")+ coord_cartesian(ylim=c(0,1.5))

Results Summary
When split up by sex, HOR/NOR RRS < 1 for males in 2 of 4 years and for females in 3 of 4 years.
HOR/F1 RRS < 1 for males in 2 of 4 years and for females in 3 of 4 years.
F1/NOR RRS < 1 for males in 1 year and zero years for females.